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We derive exact Langevin-type equations governing quasispecies dynamics. The inherent mul- 
tiplicative noise has both real and imaginary parts. The numerical simulation of the underlying 
complex stochastic partial differential equations is carried out employing the Cholesky decompo- 
sition for the noise covariance matrix. This noise produces unavoidable spatio-temporal density 
fluctuations about the mean field value. In two dimensions, the fluctuations are suppressed only 
when the diffusion time scale is much smaller than the amplification time scale for the master species. 
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Deterministic descriptions of reacting and diffusing chemical and molecular species fail to account for the system's 
internal fluctuations. Nevertheless, it is known that if the spatial dimensionality d of the system is smaller than a 
certain upper critical dimension dc, these intrinsic fluctuations can play a crucial role in the asymptotic late time 
behavior of decay rates (anomalous kinetics) and the results obtained from the mean field equations are not correct U| . 
Even far from the asymptotic regimes, these fluctuations can also control the dynamics on local spatial and temporal 
scales 0- The mean field result is only valid in the idealized infinite diffusion limit, because the reactions themselves 
induce local microscopic density fluctuations that must be taken into account in the underlying nonlinear dynamics. 

The proper inclusion of the effects of microscopic density fluctuations in reaction-diffusion systems can be carried 
out once the microscopic kinetic equations are specified. With the reaction scheme in hand, one can derive the 
corresponding continuous-time master equation, then represent this stochastic process by second-quantized bosonic 
operators and in the final step, pass to a path integral to map the system onto a continuum stochastic field theory 3, 4j. 
This technique has opened up the way for employing powerful field-theoretic renormalization group (RG) methods 
for studying fluctuations in a number of simple reaction-diffusion problems 0. Moreover, effective Langevin-type 
equations can be deduced from these field-theoretic actions, in which the noise is made manifest and is specified 
precisely. Langevin-type equations are ideally suited for investigating problems in stability and pattern formation, 
and lend themselves for the direct numerical solution of the dynamics. However, most of the studies devoted to 
fluctuations in reaction-diffusion systems are based on applying RG methods to the field-theoretic actions, with little 
attention being payed to the analytical or numerical study of the associated effective Langevin equations. This is most 
likely due to the fact that the noise in this latter representation is often imaginary or even complex, a feature, which 
at first glance, may be somewhat surprising 0, and has presented a challenging problem for numerical simulation. 

The purpose of this paper is to confront this imaginary/complex noise issue at face value. We propose and test 
out an algorithm for numerically integrating complex noise in multi-component reaction-diffusion equations |3| ■ The 
model we treat can be considered as a single quasispecies with error tail, which is analytically and computationally 
amenable. The current great revival of interest in quasispecies dynamics owes to the fact that viral population 
dynamics is known to be described by quasispecies Q • The major par t of this work has been devoted to the analysis 
of their dynamics under spatially homogeneous conditions [§. lol IllL . More recently, the importance of diffusive 
forces has been recognized and taken into account , but rather less attention has been payed to the presence 

of the unavoidable internal density fiuctuations [13 that are necessarily present in all realistic incompletely mixed 
diffusing systems of reacting agents |l6j . In view of the above considerations, it is important to understand how 
internal fluctuations affect the evolution of replicator dynamics, and in what ways do the deterministic and stochastic 
effects compete. The specific model treated here maps exactly to a set of Langevin equations. The advantage for 
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the numerical simulation is that the model has few fields, the corresponding noise terms are known exactly, and no 
approximation is required. 

In Section m we introduce the specific reaction scheme. Following a well-established procedure 0, we derive a 
field theoretic description of these reactions by means of the Doi-Pcliti formalism 01 ■ We obtain the continuum 
action, and from this derive an equivalent and exact Langevin equation description of this quasispecies model. The 
advantage of this is that the noise properties are specified automatically and indicate how the naive mean field reaction- 
diffusion equations must be modified to take into account properly the (unavoidable) internal density fluctuations. The 
ensuing noise is complex and multiplicative, and in magnitude is controlled by the competition between the replicator 
amplification and diffusion. Their numerical solution is rendered possible by employing the Cholesky decomposition 
for the associated noise covariance matrix, as we describe in detail below in Section IITTI In Sec II VI we present results 
of the numerical simulations of the complex Langevin equations derived in Sec lIIII Conclusions are drawn in SeclVI 



II. FROM THE REACTION SCHEME TO THE STOCHASTIC PDES 



A. The Model 



We consider a simple replicator model with error introduced via the faulty self-replication into a mutant species. 
The mutant species or, error-tail, undergoes non-catalyzed self-reproduction, but has no effect on the main species. 
The system is closed, only energy can be exchanged with the surroundings, where activated monomers react to build 
up self-replicative units. These energy rich monomers are regenerated from the by-product of the reactions by means 
of a recycle mechanism (driven by an external source of photons-sunlight) maintaining the system out of equilibrium. 
The closure of the system directly imposes a selection pressure on the population. In what follows. A/*, /, le denote 
the concentrations of the activated energy rich monomers, the replicators and the mutant copies, respectively. The 
kinetic constants are introduced in the reaction steps as follows: 

Accurate noncatalytic replication with rate A: 



M*+I^2I. (1) 



Error noncatalytic replication: 



M*+I^^^^^I + Ie. (2) 



Error-species replication with rate A^, 



M*+I, ^2Ie. (3) 
Species degradation and subsequent monomer reactivation with rates r, rg-. 

I^M*, (4) 
le^M*. (5) 

The quality factor Q e [0,1]. In order to keep the following development mathematically manageable, we will 
assume that the monomer reactivation step proceeds sufficiently rapidly so that we can in effect, regard the decay of 
/ and le plus the subsequent reactivation M ^ yi* ^s occurring in one single step as indicated in Eq. H4I5|I . If we 
suppose the system is being bathed continuously by an external energy source, the monomer reactivation is occurring 
continuously, and this should be a reasonable approximation. To complete the specification of the model, we will 
include spatial diffusion. We allow the M*,I and le particles to diffuse with constants Ds,Dj and De, respectively. 
Diffusion is incorporated at the outset in the master equation. The constraint of constant total particle number 
is automatically satisfied by the continuous chemical fields in the mean-field limit, as we demonstrate below. Most 
importantly, this constraint provides a selection pressure on the quasispecies. In the following, we keep the dependence 
on the model parameters general, though later on we will choose Dg >> Dj = De and r = Ve- 
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B. Mapping to bosonic field tlieory 



This chemical master equation for the model Eas. (|ll2l3l4l5l) can be mapped to a second-quantized description 
following a procedure developed by Doi 3]. Briefly, we introduce annihilation and creation operators a and 
for M* , b and for / and c and for If, at each lattice site, with the commutation relations [ai,a|] = Sij, 
[bi,b]j] — Sij and [ci,c]] = Sij. The vacuum state (corresponding to the configuration containing zero particles) 
satisfies a,i|0) = bi\0) ~ Ci\0) = 0. We then define the time-dependent state vector 



{k},{m},{n} i 



(6) 



where P{{k}, {m}, {n}, t) is the probability distribution to find fc, m, n particles of type M* , I, /g, respectively, at each 
site. The master equation can then be written as a Schrodinger-like equation 



dt 



(7) 



where the lattice hamiltonian or time-evolution operator is a fimction of a, , b, b\ c, and is given by 



aqE[« 



tbjblbi - ala^blbi 



^(1 - E h^l^ 



(8) 



Now take the continuum limit {I 0), and obtain a representation as a path integral Q over continuous fields 
a{x, t), a* {x,t), b{x, t), b* {x,t), c{x, t) and c*{x, t) with a weight exp { — S[a,a* ,b,b* ,c,c*]) , whose action S is given by 



5* 



dtd'^x 



*dta + DsWa*Wa + b*dtb + DiVb*Vb + c*dtc + D^Vc* Vc - AQ{ab*^b ~ a*ab*b) 



-A{1 — Q){ab*bc* — a*ab*b^ — Ac(^ac*'^c ~ a*ac*c) — r{a*b — 6*6) — re(a*c — c*c) 



(9) 



The stationarity conditions (the "classical field equations") SS/Sa — SS/Sb — SS/Sc — and SS/Sa* 
SS/Sc* = yield, respectively, a* = 6* = c* = 1 and the usual mean-field rate equations 

dta — DgV^a ~ Aab ~ A^ac + r6 + r^c 

dtb = DiV'^b + AQab-rb 

dtc = De'^^c + A{l-Q)ab + Aeac-reC. 



SS/Sb* 



(10) 

(11) 
(12) 



We emphasize that these equations represent the mean field approximation wherein all fluctuations are simply ignored. 
The exact and correct dynamical equations are fully stochastic, and we derive them below. Note in the mean fleld 
limit, the total particle number N is automatically conserved: adding up Eq. H10lllll2f) allows us to prove that for a 
closed reaction system (in a bounded and closed reaction domain) 



dN _ d 
'dt ^ di 



d'^x(^a{x, t) + b{x, t) + c{x, t)j = 0, 



(13) 



so that iV is a constant. 



C. Equivalent Langevin equation description 

Here, we go beyond the mean-field approximation Eq. H10I11I12|I and obtain the exact stochastic partial differential 
equations that govern the quasispecies dynamics. To do so, shift conjugate fields as follows a* = 1 + a, 6* = 1 + 6, 
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and c* = 1 + c, then we can write the action as follows: 



S = J dtd'^x [a(dta- Ds\7'^a + Aab + Aeac- rb 



b(dtb - Di^^b - AQab + rb) + d{dtc - De^^c - A{\ ~ Q)ab - A^ac 
AQabib^ - db) - A{1 - Q)ab(bc - db) - Aeac{c? - dc) 



+ rec) 



(14) 



The next step is to introduce Gaussian-distributed noise fields ria{x, t), 'q}j{x, t), r]c{x, t) which will permit us to integrate 
over the conjugate fields and obtain the exact and equivalent Langevin representation of the stochastic dynamics 
contained in S. Now the part in — S* quadratic in the conjugate fields a, 6, c contributes to the exponential weight the 
following expression (we suppress writing out the x, t dependence and the integrals J df^x dt ; these are understood to 
be included in what follows): 



exp(— 5) I quadratic = exj>\^ + AQabb — Aabdb + A{1 - Q)abbc + Aeac[c 
= exp(^ + S-V-S 
= expy^SiVj-jSj, 



;;2 



ac 



(15) 



where the vector S = (a, 6, c) and the 3x3 array V is given by 





^Aab 



'^Aab 
-AQab 



-\Aeac +\A{\ - Q)ab 
Now make use of the Hubbard-Stratanovich transformation: 



^ Aq ac 
\A{l-Q)ab 

-\-Apac 



j \{dr^, expi 



Vii^ij ^)Vj + X! "^^^i ( ^ constant x exp 



(16) 



(17) 



to express the righthand side of Eq. (|15|l as an integral over noise fields (77^). The covariance matrix V is actually 
a 3 X 3 matrix in field space and is proportional to space and time delta functions (infinite dimensional continuous 
"matrices", etc.) We immediately read off the direct and crossed noise correlations directly from Vij, since 



v= ( ivb-na) ivbVb) ivbVc 

(VcVa) iVcVb) iVcVc 



(18) 



For the final step we use Eq. H17|) to replace the right hand side of Eq. (|15|l in Eq. (|14|l . We can now integrate 

exactly over the conjugate fields d,b,c, appearing in the path integral J VaVdVb'Db'Dc'Dce'^^'^''^ '''''''^''^^ which yields 
a product of delta-functional constraints which imply the following set of exact coupled set of Langevin equations: 



dta ~ Ds\7^a ^ Aab ~ A^ac + rb + r^c + rja, 
dtb = Di'^^b + AQab~rb + r]b, 
dtc = DeS/^c + A{1 ~ Q)ab + Aettc — reC + r]c 
with noise correlations (compare (|16|l with (|18|) ) 



(19) 
(20) 
(21) 



{Va{x,t) 
{r]a{x,t)T]a{x',t') 

{r]b{x,t)r]b{x' ,t') 
{r]c{x,t)ricix',t') 
{r]a{x,t)r]b{x' ,t') 
{r]aix,t)r]c{x',t') 
{Vb{x,t)f]c{x\t') 



(r,,(x,t)) =0 



{Vb{x,t)) 


+AQa{x, t)b{x, t)5'^{x - x')S(t ~ t') 

+Aea{x, t)c{x, t)5'^{x - x')5{t - t') 

- \Aa{x, t)b{x, t)d''{x - x')S{t - t') 

-\Aea{x, t)c{x, t)5'^{x - x')S{t - t') 

+ \A{1 - Q)a{x, t)b{x, t)S'^{x - x')S{t - t'). 

Note the noise r^a has finite cross correlations Eas. l|2t)l27ll but zero autocorrelation Ea. (|23|) . This is already a 
indicator that the noise cannot be purely real. In the following Section, we will show that the pattern of the above 
correlations is solved by complex noise. 



(22) 
(23) 
(24) 
(25) 
(26) 
(27) 
(28) 



5 



III. LANGEVIN EQUATIONS WITH COMPLEX NOISE 



It is reasonable to assume that both the rephcating and mutant species diffuse with equal rates I?/e = Dj = D 
and have equal degradation rates, i.e. r^/r — 1. We thus consider the following reaction-diffusion system, in a 
two-dimensional space, subject to noise and employing the dimensionless fields, noises and model parameters (for the 
details of the non-dimensionalization of the stochastic reaction-diffusion system, see Appendix : 



da 
dh 
dc 

a7 



ab 



ac + b+ -c^ 





(29) 



V^c -I- (5(1 - Q)ab + 5ac-c + fjc 



where 



and f] = (?7f,, ?7c, ?7a) is the noise vector defined above. The initial condition ao, 6o, cq and the ratio 



of replication rates 6 = ^ < \ obeys the dimensionless constraint for the closed system N = J J dxdy{ao -I- 6o -I- ^). 
In the two-dimensional case the total number of particles is given by the ratio N — N/e where e — A/Dj is the ratio 
of the reaction to the diffusion processes (in any dimension d, we have e = (r/DiY^^A/r). 

In the deterministic case and per each single cell Sx x Sy, this reaction-diffusion system Eq. (|29|l has the following 
set of homogeneous and static solutions: 

• b ^ c — 0, a — N, absorbing solution 

• 6 = 0, c/^ = TV - i, a = i, if (5 > 1/iV 

•-c/S^^%^, b^i^^0§^,a = l/Q ifg>^andQ>l/iV. 
For convenience we can reorder the noise vector components ff = (f)b, f}c, fja) such that it has the correlation matrix: 



B =< Tfrf^ >= I ie<5(l 



eQdb _ ^eS{l~Q)ab 
Q)db 



eSac 



^eab 



'^eSac 



(30) 



where the zero autocorrelation term is located in the last column and last row. Notice that i? is a symmetric matrix 
with det-B — —^a^bce^S'^{b + c), i.e. it is negative definite. 

For an Mx M symmetric matrix B, one can apply the Cholesky decomposition B ~ LLF to extract the square root of 

the matrix in the form of a lower triangular matrix L with La = \J Bu — X]fc=i -^5; ^^"^ ^ji — 'Z~(^n ~X]a:=i LikLjk) 
and j — i + 1, Af. This decomposition is used when the symmetric matrix is positive definite. We have applied 
this algorithm to our case with negative definite correlation matrix and we obtain then the matrix "square root" L 
where some of the terms are manifestly imaginary: 



L 








2VQ 
1 (l-Q)b-2Qc 
2\/Q ^JiQc-(l-Q)H 



\ 



(31) 



We will use this decomposition to relate the noise to a new real noise ^, a white Gaussian noise with < >= 1 (i.e 



with uncorrelated real components) such that fj — (and thus rf^ = ^ L^). We thus are able to write the internal 
noise as a linear combination of white noise terms. Notice that by doing so the condition < rfrf^ >=< LS^^'^L'^ >— 
L < ^^'^ > L"^ = LL^ = B is satisfied. This transformation will allow us to separate the real and imaginary parts of 
the noise, a useful feature to have for setting up a numerical simulation of the system in Eq. H29|) . 

Given three uncorrelated (real) white Gaussian noises ^i, ^2 and ^3 (the three components of 5) we recover 77;,, fjc 
and fja with the specified cross-correlations Eqs. (|A11IA16|I by writing: 



rib = Vea \/Qb£,i 




2\/Q y/4:Qc-{l-Q)^b - 



6 



(32) 



V4Qe-(l-Q)2b 
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The noise rfa for the nutrient field always has an imaginary component and, since this noise is feeding the nutrient 
field reaction-diffusion equation, the nutrient field also has an imaginary component. Finally since all the equations 
are coupled to the nutrient a, even if the initial condition for a, b, c is real, the fields will, in principle, have both 
real and imaginary part and thus we have to solve a system of six partial differential equations, one for each field 
(Re(a), Im(a), Re(6), Im(6), Re(c), Im(c)) with two-dimensional diffusion and noise. At this juncture, it is important 
that we confirm numerically that the complex solutions, once averaged over the fluctuations, do indeed yield real 
results, in accord with the theoretical expectations [l^ . 



IV. NUMERICAL RESULTS 



By inspection of the correlation matrix we see that the internal, unavoidable, reaction noise is multiplicative 
and its intensity is proportional to e. The parameter e is the ratio of the reaction to the diffusion processes (in d = 2, 
e — A/Di and e5 = Ae/De). The multiplicative noise is proportional both to a, and either 6 or a combination of 
h with c. If the system stays close to the mean field result, then a — 1/Q and h and c scale as {NQ ~ 1). Finally, 
the multiplicative noise will increase with both increasing e and/or increasing N . Only when the diffusion processes 
dominate and e — > 0, does the noise terms vanish, the local details of reaction are erased, and we expect to recover 
the homogeneous solution of the mean-field approximation. Our main interest here is the limit when the diffusion 
does not dominate, e > and the noisy multiplicative term has a significant contribution. In this regime the problem 
can thus not be analyzed by perturbation theory and has to be treated numerically. 

The numerical simulations of system evolution have been performed using forward Euler integration of the finite- 
difference equations following discretization of space and time in the stochastic partial differential equations. The 
spatial mesh consists of a lattice of 154 x 154 cells with cell size Ax = Ay = 0.35 and periodic boundary conditions. 
Noise has been discretized as well. The system has been numerically integrated up to r = 100 (with time step 
At = 2.5 X 10""*). Integrating the system numerically we confirmed that the imaginary fields were zero in average 
as expected, since the stochastic averages (a), (&), and (c) correspond to the physical densities and thus the scaled 
number of particles was preserved and remained real: Re(iV) = N, and Im(7V) — (within computational errors) |20j| . 

This method is able to provide the time evolution of the spatial distribution of each of the fields. As an example 
we consider the evolution of the spatially uniform initial condition, d = 1/Q,b ~ ' ^/^ = ^^'cg(i-l)~'^^ 

(which is the solution of mean-field problem) for Ds/D = 10, Q = 0.92,5 = 0.8 and N — 100 and noise intensity 
e — 1. In Fig. ^ we show the time evolution of the real part of the distributions of the master (6) and mutant 
(c) species, which exhibit spatial fluctuations with respect to the mean-field value (up to +1% in black and — 1% 
in white). The spatial distributions of the master and mutant species are anti-correlated. After a transient time, 
the short-range spatial fluctuations (e.g., at r = 0.5) evolve into long-range spatial fluctuations (e.g., at r = 75) 
with wider regions of higher densities of either one or the other species. The system converges always to the same 
spatially averaged solution independently of the initial condition and noise intensity. For example, in Fig. |21-(Left) 
we show for Ds/D — 10, Q — 0.92,6 — 0.8 and N = 100, the time evolution of the spatially averaged solution 
(< a(r) >= J J dxdy a{x,y,T), etc.) for two noise intensities, e = 0.1 and e = 10~^, when the initial condition is set 
to a (a;, ij, 0) — b{x, y, 0) = c{x, y,0)/5 = N/3. We can see that independently of the noise intensity this averaged value 
converges to the solution of the mcan-fleld problem d = 1/Q,b — c/J — ^^'0(1-^)''^'' ■ Higher order 

moments, such as the mean square value, can be also extracted from the distribution. For the same study case we 
the show the time evolution of the mean square value (cr^^^ = a/ < (6(t)— < 6(t) >)2 >) for three noise intensities, 
e — 1, e — 0.1 and e — 0.01, Fig. [21-(Right), conflrming that only in the limit e ^ do spatio-temporal fluctuations 
vanish and we recover the homogeneous stationary solution of the mean-field approximation. 



V. DISCUSSION 



Starting from the microscopic kinetic equations for a single master species competing with mutants for a limited 
pool of nutrients, we have derived a set of stochastic partial differential equations that exactly describe the complete 
dynamics of this closed reacting and diffusing system. The unavoidable fluctuations inherent to the system give rise 
to multiplicative noise having both real and imaginary components. A procedure is presented to solve numerically 
this set of partial differential equations. The internal noise associated with the microscopic details of the reaction 
produces unavoidable spatio-temporal density fluctuations around the mean field value. We estimate the size of these 
fluctuations. They strictly vanish, and the mean-field limit is recovered, only when the diffusion processes are much 
faster than the rate of master species amplification. We find that the averages {d),{b) and (c) tend to stationary 
values (see Fig|2KLeft) ) as do also the second moments, or variances, aa{T):'^b{T) ^^'^ '^cir), see Figl^KRight). 
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FIG. 1: Temporal sequence of spatial fluctuations in the real part of the stochastic fields b{x,y,T) (upper row) and c(f,y,r) 
(lower row). The spatial distributions are shown in (a) at r = 0.5, in (b) at r = 50 and in (c) at r = 75. 




FIG. 2: (Left) Time evolution of spatially averaged values (upper curve for < 6 >, middle for < c > /5 and lower for < a >) 
for A'^ = 100, e = 0.1 and e — 0.01 (the difference in these two noise levels is not distinguishable in the figure) when the 
initial condition was set to a{x, y, 0) — b{x, y, 0) = c{x, y, 0)/S = N /S . The spatially averaged solution always converges to the 
mean-field solution which is marked by the horizontal lines. (Right) Time evolution of the mean square value crj,^^) for e — 1, 
e = 0.1 and e = 0.01 showing that the spatial fluctuations scale with the noise intensity, vanishing only in the limit e — > 0. 



The general purpose algorithm presented here expresses a set of complex Gaussian noises with defined covariance 
matrix as a linear combination of real, white Gaussian noises. This allows for the numerical generation of this noise, 
and thus the numerical integration of Langevin-type equations with complex noise. We have performed numerical 
simulations with these dynamical equations to assess directly the influence of this noise on the evolution of the 
stochastic fields associated with a quasispecies, its error tail and the activated monomer distributions. Apart from 
the specific application to the problem analyzed in this paper, we believe this matrix decomposition will be of great 
practical use for simulating other stochastic PDEs with complex noise, a subject which has received little attention 
up to now 7J. In this regard it is interesting to emphasize that the appearance and necessity of complex and/or 
imaginary noise in probabilistic descriptions of both quantum optics and nonlinear chemical reaction systems has 
been recognized for some time now and began to be placed on a firm theoretical footing well over twenty years ago 

The simulations are carried out in d = 2 space dimensions. For this dimension, the basic amplitude of the noise 
is controlled by the dimensionless parameter e = A/Dj, where A denotes the amplification rate of the quasispecies 
and Di its diffusion. Thus, the noise is controlled, in part, by the competition between production and diffusion. 
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Diffusion tends to erase or smooth out local concentration gradients while particle production, which occurs locally, 
tends to increase them (i.e., increases the local density of the species being amplified). In d = 2 and using Eq. HA1|I . 
we can also write e = t^/tA, where 1^,1 a denote the diffusion and species amplification time scales, respectively. This 
shows that the noise arises through the competition between these two time scales. Furthermore, since the noise is 
multiplicative^ its amplitude also depends on the bilinear product of concentration fields (see, e.g., Eqs. I|A10I - IXT6|) ) 
so that whenever monomer and replicator or monomer and error tail meet at a point and react, this gives rise to noise 
at that point, whose strength is directly proportional to the product of particle concentrations at that point (equal to 
the product of particle numbers per unit area). Thus the complete noise amplitude is modulated by e and the local 
bilinear concentrations. So, while diffusion smooths out inhomogeneities and tends to diminish the internal noise, 
increasing the total fixed particle number N (in a bounded domain) leads to stronger local fluctuations. This is easy 
to understand since, in spite of the diffusion, increasing N leads to a "pile-up" of reactants at spatial points. There 
will always be fluctuations about this approximate homogeneous state, and the statistical deviation from homogeneity 
grows with increasing total particle number. 

In realistic situations with finite diffusion rates, the system does not converge to a homogeneous solution with a 
unique defined value, but to a state where the densities fluctuate both in space and time around a mean value. This 
approach is useful to estimate the expected deviations of the densities with respect to the mean-field value when the 
microscopic reaction details are taken into account. For systems with higher degree of non-linearity these deviations 
may eventually lead the system to new asymptotic states and also induce the formation of true patterns. We point out 
that no true spatial patterns are generated by the underlying reaction diffusion model studied in this paper. The model 
as it stands is weakly non-linear, depending only quadratically on the fields and these second-order reaction terms are 
not capable of giving rise to spatial patterns. The proper inclusion of noise induces random structures. However, once 
catalyzed self-replication is included, thereby leading to a bona-fide network of quasispecies, the underlying dynamics 
becomes cubically nonlinear, and cubic terms can lead to a variety of possible spatial patterns whose evolution 
and stability properties in presence of internal noise can be profitably studied with the methods offered in this paper 
[THj l. Finally, we mention that the methods in were used recently to study the critical behavior of a simple 

model of quasispecies in the vicinity of the error-threshold |l9j . 
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APPENDIX A: NONDIMENSIONALIZATION 



For the purposes of numerical simulation, it is convenient to cast the system of stochastic partial differential 
equations Eqs. H19I2I|I and noise statistics Eqs. (|22l28|) in terms of dimensionless fields and parameters. To this end, 
define basic length and time scales L and T, respectively. Then dimensional analysis of any of the three stochastic 
equations in Eqs. (|I9I21|) yields 

[a]^[h] = [c\=L-\ [Di] = [D,] = [Ds\=L^/T [Q] = I, [A] = [A,] ^ L'' /T, (Af) 
[r] = [r,] = T-\ [i^,] = [r?,] = [r?,] = T'^L-^. (A2) 

Define the dimensionless fields: 

a= — a, h — —h, c= — -c, (A3) 
r r r 

and the dimensionless time and spatial coordinates and their corresponding derivative operators: 

/ r 

T = rt, xj = i^ — j Xj, (A4) 

did - , /f/x^o , . ^, 

^ 7r = -^' V^ = {—)V^. A5 
OT r at \ J. ' 
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We obtain the following dimensionless version of the stochastic equations listed in Eqs. (|19I21|I : 



ab - 



drb = V^ + Qdb-b + ijb, 



my 



\DiJ 

where the dimensionless noises are defined by 

. _ A 

and the dimensionless noise correlations are 



A 



[l~Q)db- 



Ae 

A 



ac - 



A 

—Vb, 





iVaix 


,r)) = 


{Vb{x,T)) 


= {flc{x,T)) = 









T)fla{x , 


r')) = 











iVbix 


T)flb{x\ 


r')) = 


+e Qa{x, 


T)b{x,T)S'^{x- 


-x')SiT- 


-r') 


{Vc{x 


T)flc{x\ 


r')) = 


+6(5^ a{x, 


t)c{x,t)S'^{x - 


-x')5{t^ 






T)flb{x\ 


r')) = 


— iea(x, 


T)b{x,T)S'^{x - 


-x')SiT- 


-r') 


{Va{x 


T)f]c{x\ 


r')) = 


— ^£(5 d{x 


, t)c(x, t)S'^{x 


-x')S{t 


~r') 


{Vb{x 


T)f]c{x', 


r')> = 


+^eS{l- 


- Q)d{x, T)b{x, 


T)6'^ix- 


x')S{ 



In arriving at this, note that 



[S'{x)] = L- 



whereas the following deltas are dimensionless: 

[S'ix)] 

and 



= 1, 



m] = T- 



6''ix)6{t) = r 



d/2 



The basic control parameters are given by 

A( r 



\DjJ 



6'^{x)6{t). 



< 1, 1 > Q > 0. 



(A6) 
(A7) 
(A8) 

(A9) 



(AlO) 
(All) 
(A12) 
(A13) 
(A14) 
(A15) 
(A16) 



(A17) 



(A18) 



(A19) 



(A20) 



Note that in d = 2 space dimensions the noise amplitude e is determined uniquely by the ratio of the replication rate 
A to replicator diffusion Dj. It is reasonable to assume that both master sequence and error tail replicator molecules 
decay with the same rate, r — r^. Moreover, since both the master sequence and error tail are built up from the same 
monomer pool and have the same total length, they should also have identical diffusion constants Dj = Dg. The 
much smaller monomers should diffuse more rapidly, so we can take Dg >> Dj ^ = D . Lastly, we note that the 
noise-averaged fields satisfy the constraint: 



Syi{d + b+\c) ^ N ^ tN, 



where N is the total particle number. 

The complete nondimensional model has five parameters: e,6,Q, D^/ D > 1 and N. 



(A21) 



[1] K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 (1984). 

[2] M.R Zorzano, D. Hochberg and F. Moran, Physica A 334, 67 (2004). 

[3] M. Doi, J. Phys. A: Math. Gen. 9:1465, 1479 (1976). 



10 

[4] L. Peliti, J. Physique 46:1469 (1985). 

[5] For a recent review, see U.C. Tauber, M. Howard and B.P. Volniayr-Lee, J. Phys. A: Math. Gen. 38, R79 (2005). 
[6] M.J. Howard and U.C. Tauber, ,]. Phys. A:Math. Gen. 30, 7721 (1997). 

[7] A single imaginary noise Langevin equation without diffusion has been treated in O. Deloubriere, L. Frachebourg, H.J. 

Hilhorst and K. Kitahara, Physica A 308, 135 (2002). 
[8] E. Domingo, O.K. Biebricher, M. Eigen and J.J. Holland, Quasispecies and RNA Viruses: Principles and Consequences 

(Eureka.com, Austin, 2001). 

[9] M. Eigen and P. Schuster, The Hypercycle-A Principle of Natural Self-Organization (Springer,Berlin, 1979). 
[10] C. Biebricher and M. Eigen, Virus Research 107, 117, (2005). 

[11] C.K. Biebricher and M. Eigen, Current Topics in Microbiol. Immunol. 299, 1 (2006). 

[12] M.A. Andrade, J.C. Nufio, F. Moran, F. Montero and G.J. Mpitsos, Physica D 63, 21 (1993). 

[13] P. Chacon and J.C. Nuno, Physica D 81, 398 (1995). 

[14] M.B. Cronhjort and C. Blomberg, Physica D 101, 289 (1997). 

[15] D. Hochberg, M.-P. Zorzano and F. Moran, J. Chem. Phys. 122, 214701 (2005). 

[16] I.R. Epstein, Nature 374, 321 (1995). 

[17] J.L. Cardy, in The Mathematical Beauty of Physics (World-Scientific, Singapore, 1996) eds. J.-M. Drouffe and J.-B. Zuber. 
[18] P.D. Drummond, C.W. Gardiner and D.F. Walls, Phys. Rev. A24, 914 (1981). 
[19] R. Pastor-Satorras and R.V. Sole, Phys. Rev. E 64, 051909 (2001). 

[20] In fact, for computational purposes, Im(A'') was set as small as possible. We confirmed that both Re(Af) = / / d£-dy(Re(a) + 
Re(6) + 5£^) ^ ff and Im(iV) = / / didy (Im(a) +Im(6) + l51^) = N x 10"'^ are conserved during the integrating time 
and that the results were independent of Im(A''). 



